#!/usr/bin/perl -w
use strict;
use FindBin;
use lib ("$FindBin::Bin/..", "/net/cpp-group/Leo/bin");
use cluster_from_pairs;
use parse_bl2seq;
use DBI;
use DBD::Pg;
use linewrap;
use Getopt::Long;

#88888888888888888888888888888888888888888888888888888888888888888888888888888888

#	Usage

#88888888888888888888888888888888888888888888888888888888888888888888888888888888

my $usage = <<'USAGE';

USAGE:

    write_species_sequences_identifiers.pl
                --species e_g_canis
                --directory PATH
				--help
	
	Description:
		Writes sequences and identifier map for 
		species pair to a specied directory
USAGE


#88888888888888888888888888888888888888888888888888888888888888888888888888888888888888888

#	Get options

#88888888888888888888888888888888888888888888888888888888888888888888888888888888888888888
# mandatory

# optional parameters
my $help;
my $species;
my $directory;
GetOptions(
			'species=s'				=> \$species,
			'directory=s'			=> \$directory,
			'help'					=> \$help);

die $usage if ($help);
die $usage unless ($species && $directory);


#-----------------------------------------------------------------------------------------
#	Get DB parameters
#
use db_parameters;
#
#	Get DB parameters
#-----------------------------------------------------------------------------------------






#-----------------------------------------------------------------------------------------
#
#	Get prot sequences
#
#
sub get_all_prot($)
{
	my ($dbh) = @_;
	my $sql_cmd =<<"PL/SQLCMD";
    SELECT
        DISTINCT ON (cdna_hash)
        cdna_hash, sequence
    FROM
        cdna_sequences 
        JOIN sequences USING (pi) 
        JOIN ens_id USING (prot_id)
    WHERE
        species = '$species'
PL/SQLCMD
	return $dbh->selectall_arrayref($sql_cmd)		or die DBI::errstr;
}



#-----------------------------------------------------------------------------------------
#
#	Get cDNA sequences
#
#
sub get_all_cdna($)
{
	my ($dbh) = @_;
	my $sql_cmd =<<"PL/SQLCMD";
    SELECT
        DISTINCT ON (cdna_hash)
        cdna_hash, cdna
    FROM
        cdna_sequences 
        JOIN ens_id USING (prot_id)
    WHERE
        species = '$species'
PL/SQLCMD
	return $dbh->selectall_arrayref($sql_cmd)		or die DBI::errstr;
}

#-----------------------------------------------------------------------------------------
#
#	get ensp to ensg
#
#
sub get_taxon_identifier_map($)
{
	my ($dbh) = @_;
	my $sql_cmd =<<"PL/SQLCMD";
	SELECT
        DISTINCT ON (prot_id)
		prot_id, gene_id, species
	FROM
        ens_id
    WHERE
        species = '$species'
PL/SQLCMD
	return $dbh->selectall_arrayref($sql_cmd)		or die DBI::errstr;
}

sub get_hash_to_prot_id_map($)
{
	my ($dbh) = @_;
	my $sql_cmd =<<"PL/SQLCMD";
	SELECT
        DISTINCT ON (prot_id)
		cdna_hash, prot_id
	FROM
        cdna_sequences 
        JOIN ens_id USING (prot_id)
    WHERE
        species = '$species'
PL/SQLCMD
	return $dbh->selectall_arrayref($sql_cmd)		or die DBI::errstr;
}


sub get_hash_to_gene_id_map($)
{
	my ($dbh) = @_;
	my $sql_cmd =<<"PL/SQLCMD";
	SELECT
        DISTINCT ON (cdna_hash)
		cdna_hash, gene_id, species
	FROM
        cdna_sequences 
        JOIN ens_id USING (prot_id)
    WHERE
        species = '$species'
PL/SQLCMD
	return $dbh->selectall_arrayref($sql_cmd)		or die DBI::errstr;
}


#-----------------------------------------------------------------------------------------
#
#	get identifiers
#
#
sub get_all_cdna_hash($)
{
	my ($dbh) = @_;
	my $sql_cmd =<<"PL/SQLCMD";
	SELECT
        DISTINCT ON (cdna_hash)
		cdna_hash
	FROM
        cdna_sequences 
        JOIN ens_id USING (prot_id)
    WHERE
        species = '$species'
PL/SQLCMD
	return $dbh->selectall_arrayref($sql_cmd)		or die DBI::errstr;
}




#-----------------------------------------------------------------------------------------
#
#	check_progress
#
#
sub check_progress
{
	if (++$_[0] % $_[1] == 0)
	{
		print STDERR "\n\t" if (($_[0] / $_[1]) % 50 == 0);
		print STDERR ".";
	}

}


sub print_sequences($$&$)
{
	my ($sequence_type, $file_name, $func, $dbh) = @_;
	print STDERR "\tQuerying panda for $sequence_type $species sequences ...\n";
	open SEQ, ">$directory/$file_name" or
		die "Error\n\tCould not open $directory/$file_name\n$!";
	my $sequences = $func->($dbh);
	my $cnt_seq = 0;
	print STDERR "\t[";
	for my $seq(@$sequences)
	{
		check_progress($cnt_seq, 2000);
		print SEQ ">$seq->[0]\n", linewrap($seq->[1], 60);
	}
	print STDERR "]\n";
	print STDERR "\tFreeing memory...\n";
}

sub print_identifiers($&$)
{
	my ($file_name, $func, $dbh) = @_;
	print STDERR "\tPrinting $species identifiers to $file_name ...\n";

	open IDS, ">$directory/$file_name" or
		die "Error\n\tCould not open $directory/$file_name\n$!";
	print STDERR "\tQuerying panda for IDs...\n";
	my $ids = $func->($dbh);
	for my $id(@$ids)
	{
		print IDS join ("\t", @$id), "\n";
	}
	print STDERR "\tDone.\n";
}



















#88888888888888888888888888888888888888888888888888888888888888888888888888888888888888888
#
#	main logic
# 



#
# connect to genome database
#
my $dbh = connect_to_panda();

#
#	print all prot & cDNA sequences
#
print_sequences("protein",  "prot.mapped.fa", \&get_all_prot, $dbh);
print_sequences("cDNA",  	"cdna.mapped.fa", \&get_all_cdna, $dbh);

#
#	print identifiers
#
print_identifiers("prot_to_gene_ids.map",				\&get_taxon_identifier_map,	$dbh);
print_identifiers("cdna_hash_to_prot_id.map",			\&get_hash_to_prot_id_map,	$dbh);
print_identifiers("cdna_hash_to_gene_id_species.map",	\&get_hash_to_gene_id_map,	$dbh);
print_identifiers("cdna.ids",							\&get_all_cdna_hash,		$dbh);


$dbh->disconnect();
